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We present a high order perturbation approach to quan- 
titatively calculate spectral densities in three distinct steps 
starting from the model Hamiltonian and the observables of 
interest. The approach is based on the perturbative continu- 
ous unitary transformation introduced previously. It is con- 
ceived to work particularly well in models allowing a clear 
identification of the elementary excitations above the ground 
state. These are then viewed as quasi-particles above the vac- 
uum. The article focuses on the technical aspects and includes 
a discussion of series extrapolation schemes. The strength of 
the method is demonstrated for 5=1/2 two- leg Heisenberg 
ladders, for which results are presented. 

PACS numbers: 75.40. Gb, 75.50.Ee, 75.10.Jm 



I. A. INTRODUCTION 

Spectroscopic measurements provide important in- 
sights in the microscopic structure of solids. Generic 
spectroscopic data contains information about energy 
bands, associated density of states, matrix elements and 
selection rules. From the theoretical point of view this 
information is indispensable in the process of formulat- 
ing and testing appropriate microscopic models. How- 
ever, a quantitative comparison of the theoretical and 
the experimental data strongly depends on the power 
of the method used to calculate the properties of the 
microscopic model. In this context high order pertur- 
bation theory has proven to be a versatile and flexible 
tool. Especially in the field of spin models a variety of 
perturbation techniques is in use (for an overview see 
Ref . [1] ) . Most of them concentrate on the calculation of 
one-particle energies and in some circumstances on the 
associated spectral weights [2-4]. However, so far the 
important information contained in the line shapes of 
spectroscopic data, i.e., the model's spectral densities as- 
sociated with the experimental observables, has not been 
exploited. The need for a quantitative calculational tool 
closing this gap is apparent. 

In the last couple of years there has been considerable 
progress in the field of high order perturbation theory. 
For a long time the methods were restricted to ground 
state energies and one-particle energies. Just recently 
we were able to quantitatively calculate two-particle ex- 
citation energies (bound states) in various spin mod- 
els [5,6]. These calculations were based on the pertur- 
bative continuous unitary transformation (CUT) method 



introduced previously [7]. This technique was used suc- 
cessfully for low-energy calculations in various spin mod- 
els before [7-9]. The key point is to construct the trans- 
formation such that the resulting effective Hamiltonian 
is block-diagonal with respect to the number of parti- 
cles. The linked cluster series expansion, an established 
high order perturbation method, has been shown lately 
to be also well suited for calculating two-particle ener- 
gies [10,11]. The basic idea is again the application of an 
orthogonal transformation which is designed to achieve a 
block-diagonal effective Hamiltonian. 

In Ref. [12] we presented an extension of the CUT 
method allowing a systematic high order perturbation 
theory for observables. In the present article we use the 
results of Ref. [12] to calculate spectral densities of multi- 
particle excitations quantitatively. The method allows to 
obtain the complete spectral information for experimen- 
tal relevant observables without any finite size restric- 
tions. The results are exact in the sense of the thermo- 
dynamic limit, i.e., each order can be calculated for the 
thermodynamic limit. By truncating the series expan- 
sion at a (high) maximum order we restrict to dynamic 
processes for which the involved particles interact within 
a certain finite distance to each other. Obtaining higher 
orders amounts to allowing larger distances. Hence, the 
scheme can be expected to work particularly well in sys- 
tems with short correlation lengths. 

Results for various spin systems have been reported in 
earlier publications (Refs. [13-21]). The article on hand 
gives the technical details necessary to apply the method. 
We like to mention that the linked cluster method men- 
tioned above was recently extended to allow for calcu- 
lating spectral densities, too [22]. It thus constitutes an 
alternative approach. 

To illustrate our approach we consider the S" = 1/2 an- 
tiferromagnetic two-leg Heisenberg ladder as an interest- 
ing and comprehensive testing ground. The Hamiltonian 
reads 

H{J^,J^) = J^H^ + J\\H\\ (1) 

= [^J_Si_iS2,i + J^ll (Si,iSi,i+i + S2,iS2,i+l)] , 

i 

where i denotes the rung and 1, 2 the leg. 

In the broad field of spin liquid systems there has been 
an ongoing theoretical interest in the spin ladder and its 
extended versions [23-37,11,38]. The model is realized in 
a number of substances [39] and there is a large amount of 
experimental data available, see e.g. Refs. [40-47]. Addi- 
tionally, the experimental evidence for superconductivity 
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in Sro.4Ca13.6Cu24O41.84 under pressure [48] has intensi- 
fied the interest. 



I. B. ARTICLE OUTLINE 

The basic concept of our approach to spectral densities 
is as follows. For a given observable O the T = spectral 
density is calculated from 

S{lo) = --lmg{cu) , (2) 

TT 

where ^(w) is the retarded zero temperature Green func- 
tion 

(3) 

Here ipo is the ground state and Eq is the ground state 
energy of the system. Since the expectation values of 
quantum mechanical obscrvables do not change under 
unitary transformations, the Green function Q, and thus 
S, will not be altered if the operators H and O and the 
state IV'o) appearing in g are substituted by the effective 
operators (state) Hes and Oeff (|V'o,eff)) obtained from 
the CUT method. Our procedure can be divided into 
three steps: 

1. Use the CUT to derive an effective Hamiltonian 
Hes unitarily linked to H. 

2. Use the same transformation to derive the effective 
observable O^s from some initial observable O of 
interest. 

3. Evaluate Eq. (3) for the effective operators in terms 
of a continued fraction. 

Now, the key ingredient of our approach is to identify 
suitable (quasi-)particles which can be used to describe 
the T = physics of the system under study. We then 
proceed and construct the perturbative unitary transfor- 
mation such that the resulting effective Hamiltonian iJeff 
conserves the number of these particles. Let Q be the op- 
erator which counts the number of particles. Then the 
conservation of the number of quasi-particles reads 

[i^cff,0] = . (4) 

In Sect. H we describe in detail how this idea can be put 
to use by illustrating the procedure for the spin ladder 
example. 

Applying the same transformation to observables O 
different from the Hamiltonian leads to effective observ- 
ables Oeff' not conserving the number of particles in gen- 
eral. Their action on the ground state, as needed in the 
evaluation of the Green function (3), is characterized by 



the number of particles they inject, i.e., by the number 
of elementary excitations they excite. We thus decom- 
pose Oeff into operators injecting none, one, two and so 
on particles in the system. In Sect. Ill we again use the 
ladder example to illustrate the practical realization of 
this concept for experimentally relevant obscrvables. 

Once the effective Hamiltonian and the observables 
are obtained they are inserted into the Green function. 
Sect. IV features a detailed description of how the result- 
ing expression is manipulated to extract the correspond- 
ing (energy- and momentum-resolved) spectral densities. 
Again, the one-, two- and more-particle contributions to 
the total spectral density can be treated separately lead- 
ing to a simple and comprehensive physical picture in the 
end. 

In Sect. V we address the problem of series extrapola- 
tion, an inevitable difficulty in perturbative approaches. 
Following our approach, a very large number of quanti- 
ties has to be extrapolated simultaneously. This poses 
a difficulty which cannot always be tackled by standard 
techniques. We introduce a robust extrapolation scheme 
based on optimized perturbation theory [49] and show 
how it can be applied to extend the range of the pertur- 
bative results. 

The article is summarized in Sect. VI. 



II. TRANSFORMATION OF THE 
HAMILTONIAN 

We start by briefly explaining how the effective Hamil- 
tonian Hes is constructed from the initial ladder Hamil- 
tonian (1). (Other examples for this procedure can be 
found in Refs. [7,9] for instance.) In the subsequent sec- 
tions II. A through II. C we illustrate in detail how the 
zero-, one- and two-particle energies are calculated from 

Let us assume that the initially given Hamiltonian H 
can be formulated as perturbative problem 

H = U + xV . (5) 

In its present formulation (extensions are possible) the 
perturbative CUT method relies on two prerequisites 
calling for a band-diagonal problem as starting point: 

(A) The unperturbed Hamiltonian U must have an 
equidistant spectrum bounded from below. The 
difference between two successive levels is called an 
energy quantum or (quasi-)particle and we identify 
Q = U. 

(B) There is a number N B N > such that the per- 
turbing Hamiltonian V can be written as ^ = 
X]^=-Ar ^« where T„ increments (or decrements, 
if n < 0) the number of energy quanta by n: 
[Q,T„] =nT„. 
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We now show that the initial ladder Hamiltonian (1) 
meets these requirements. We reformulate the ladder 
problem according to 

^=H^+xH\\, (6) 

with X = J\\/J± as perturbation parameter, H± = 
H{1,0) and H\\ = H{0,1). We assume J± to be an- 
tifcrromagnctic and set J± = 1 henceforth. The limit 
of isolated rungs is the limit for which our perturbative 
treatment is controlled. 

The ground state of the unperturbed part H±^ is the 
product state with singlets on all rungs. A first excited 
state is a single rung excited to a triplet. There are 3L/2 
such elementary triplet excitations if L is the number of 
spins. The energetically next higher state is given by two 
rung-triplets and so on. The operator H± simply counts 
the number of rung-triplets and it is easily verified that 
condition (A) is fulfilled. 

For the rest of this article we identify Q = H±, i.e., the 
elementary excitations of the unperturbed part (rung- 
triplets) serve as (quasi-)particles in our treatment of the 
ladder system. Since we prefer the particle picture we call 
these elementary excitations triplons [18]. They are to- 
tal spin S* = 1 excitations appearing in three different 
variants (S'^-components -1,0,1) in contrast to magnons 
which have two variants only (S^ ± 1) in a phase of bro- 
ken spin symmetry. Both, triplons and magnons, must 
be distinguished from spinous, which are S = 1/2 excita- 
tions. Whenever we refer to the ladder system we will use 
the term triplon. In more general discussions we retreat 
to the term (quasi-)particle. 

As soon as we turn on the inter-rung interaction {x > 
0) the triplons become dressed particles. The central idea 
of the CUT approach is to map the initial problem onto 
an effective Hamiltonian for which the simple triplon- 
states, originally defined for the unperturbed part, can 
be used to calculate all energy levels of the system. 

We proceed and analyze the action of the perturbing 
part on the triplon-states. Let |n) denote a state 
with n rungs excited to triplets (n-triplon state), i.e., 
H_i_\n) = n\n). Then 

= T_2 + To + T2 , with (7) 

Ti\n) ~ \n + i) and 

E^o.±2M > (8) 

where v denotes pairs of adjacent rungs. The index u 
can also be viewed to count the bonds connecting adja- 
cent rungs. The action of the local operators %,±2{i^) on 
neighbouring rungs is given in Table I. Condition (B) is 
fulfilled with N = 2. There appear no T±i in 

The perturbative CUT is engendered by introducing 
an auxiliary variable £ G [0, 00]. The CUT gives rise to 
the flow equation (details in Ref. [7]) 

^E^ = [rj{x;i),H{x;£)], (9) 
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TABLE L Action of the operators % as defined by Eq. (8) 
on product states of adjacent rungs. Singlets are denoted 
by s and triplons by where the superscript indicates the 
magnetic quantum number. The remaining matrix elements 
can be found by using = T_n- 



which controls the flow of the Hamiltonian in the trans- 
formation process. We fix H{x; 0) = H{x) and define 
ileff(x) := H{x,oo). 

As shown in Ref. [7] the best choice for the infinitesimal 
unitary generator is (sgn(O) = 0) 

Vi,j{x;£) = sgn(Qi - Qj)Hij{x;i) , (10) 

where the matrix elements rji^j and Hi j are given in the 
eigen basis {\n)} of Q — H±. In the limit £ ^ 00 gen- 
erator (10) eliminates all parts of H{x;£) changing the 
number of particles, i.e., [H(.s,H±] — 0, and keeps the 
flowing Hamiltonian (intermediate £) band-diagonal [7]. 
The vanishing commutator expresses the fact that Heg is 
block-diagonal with respect to the number of particles. 

A perturbative realization of the transformation yields 
the effective Hamiltonian as operator series expansion 

00 

H^s{x) =H^ + J2x'' Yl C!{rri)T{rn) . 

k=l |m|=fc,M(m)=0 (11) 

Here m is a vector of dimension \m\ = k oi which the 
components are elements of {±A^, ±(7V — 1), ... ± 1, 0}. 
In the ladder case we have N = 2 and Ti = T_i = 
(cf. Eq.(7)). The operator products T{m) are defined by 
T(m) = T„ijT„i2 • • -Trn^, with Tm, as given in Eq. (8); 
k is the order of the process and M(m) := X^m^ = 
signifies that the sum of the indices vanishes which re- 
flects the conservation of the number of particles. Thus 
the action of -ffeff can be viewed as a weighted sum of 
virtual excitation processes T{m) in each of which the 
particle number is conserved. The coefficients C(m) can 
be calculated as fractions of integers (in the ladder case 
up to order k = 15). The effective Hamiltonian is thus 
an exact series expansion up to some maximum order. 

We want to emphasize that the effective Hamiltonian 
Hcs with known coefficients C(m) can be used straight- 
forwardly in all perturbative problems that meet condi- 
tions (A) and (B). The coefficients C(m) will be made 
available electronically on our web pages [50] . 

The action of the effective Hamiltonian (11) on the 
states of interest is calculated on a computer. In the 
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following subsections we illustrate how we obtain pertur- 
bative results for the ground state energy, the one-triplon 
and the two-triplon energies for the spin ladder from Hcs- 
On general grounds we showed previously [12] that H^s 
decomposes into a sum of irreducible n-particle operators 



ground state energy per site eg can be calculated in the 
thermodynamic limit on a finite minimum cluster by 



off 



(12) 



n=0 



For the problem on hand the operator H„ measures n- 
triplon energies no matter how many triplons are present 
as long as there are at least n triplons. On states con- 
taining less than n triplons the action of H„ is zero. The 
matrix elements of Hn are extensive quantities. By ex- 
ploiting the linked cluster theorem they can therefore be 
calculated perturbatively for the infinite system on finite 
minimum clusters, which are just large enough to per- 
form the calculations without finite size effects. More 
details can be found in Ref. [12]. 

The Hn are calculated recursively from H^s starting 
with Ho (Eqs. (9) in Ref. [12]). In this way, Eq. (12) 
stands for a systematic energy-calculation scheme. One 
starts by calculating the ground state energy (i?o) and 
proceeds by calculating one-triplon energies (Hi). True 
two-triplon interactions can be calculated by including 
i?2 and so on. A detailed description of this issue, in 
particular how the iJ„ are defined, is given in Ref. [12]. 

In the following three subsections we address the cal- 
culation of Hq, Hi and H2 for the spin ladder system 
separately. All necessary computational details are pre- 
sented. Particular attention is paid to the choice of minu- 
mum clusters for the ladder system. The aim is to offer 
a worked example for the interested reader. 



A. Zero Triplon: Hq 

Let |0) denote the triplon vacuum. This is the state 
where all rungs are occupied by singlets. Clearly, |0) is 
the ground state of H{x = 0) = H±. The one-triplon gap 
separates the corresponding ground state energy from the 
first excited level. In Ref. [51] we showed on general 
grounds that the particle vacuum |0) remains the ground 
state of Hcff for finite x unless a phase transition occurs 
(e.g. a mode softening at some critical value Xc). For the 
ladder system in particular, one observes that the one- 
triplon gap decreases on increasing x but stays finite for 
all < a; < 00 [52-54]. There are no phase transitions 
in this range and |0) remains the ground state. Since 
Hcff conserves the number of triplons we conclude that 
(0|iJcff(0 < X < oo)\0) is the ground state energy. The 
point X = cxD is a singular point at which the two legs 
of the ladder decompose into two decoupled gapless spin 
1/2 chains. 

Since the action of Hq on |0) coincides with the action 
of Hcff on this state (see Ref. [12]), every order of the 



eo - {0\Hcff\0)/{2N) 



(13) 



where N is the number of rungs used in the minimum 
cluster. 

We now specify the minimum cluster. At first, it is 
clear that we need a closed ladder segment. This en- 
sures that there are no end rungs, which are linked to 
the cluster by one inter-rung bond only. They would not 
contribute the same amount of energy as the fully linked 
rungs in the middle of the cluster. Fig. 1 shows a cluster 
of the ladder system which has been closed to a ring. 




FIG. 1. A periodically closed cluster of ten rungs. HcB to 
third order connects at maximum four neighbouring rungs by 
activated bonds (see text). The connected rungs are printed 
in grey. 



A close inspection of Eqs. (11) and (8) and Tab. I shows 
that Hcff connects a maximum of / -f 1 rungs on a fi- 
nite cluster of N rungs in l^^ order. In other words: A 
maximum of I bonds between neighbouring rungs can be 
activated in l^^ order. A bond 1/ is said to be activated, 
if a part of Hcff, i.e., the specific local operator 7^(i/) in 
r„ ~ of ^eff (see Eq. (8)), has acted on the 

two rungs connected by v. 

The linked cluster theorem states that only those pro- 
cesses induced by the T(m) of i?cff contribute to the 
ground state energy (and all other extensive quantities), 
in which all activated bonds are linked. Processes involv- 
ing disconnected active-bond distributions cannot con- 
tribute. The basic argument is sketched in Fig. 2. This 
means in our case, that a cluster of / -I- 1 rungs is suf- 
ficient to calculate the l^^ order contribution avoiding 
wrap-arounds as indicated in Fig. 3. 

Minimum number of rungs to calculate 
the l^^ order contribution to eq in the 
thermodynamic limit 

= ; + i. (14) 

Once the minimum cluster is specified it is straight for- 
ward to calculate eg. The action of i?off on |0), which 
we have calculated up to 14*'^ order, on a closed cluster 
containing 15 rungs is implemented on a computer. De- 
tails of this procedure can be found in Ref. [7] . Since Hcff 
conserves the number of triplons we have i/cff |0) ~ |0). 
The constant of proportionality is the ground state en- 
ergy Eq of the 15-rung cluster. The ground state energy 
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FIG. 2. A closed ladder-segment of six rungs. Rungs are 
depicted by circles and (active) bonds between rungs by 
(thick) solid lines. In a process of order / = 3 a maximum 
of 3 bonds can be active. On a closed cluster of AT = 6 there 
are 6 possibilities to arrange linked bonds (top row). One 
clearly sees that this number grows linearly in N. The given 
example of 3 disconnected active bonds (bottom row) has 12 
possibilities, which would lead to a super-extensive contribu- 
tion oc N'^ to the extensive quantity under study. Thus they 
do not contribute. 



per spin is finally given by eo = i?o/30. The result is a 
14th order polynomial in x. It is the exact energy of the 
infinite system to the given order 
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(15) 



The coefficients are fractions of integers and therefore 
free from rounding errors. Our findings agree with the 
numerical results given by Zheng et al. [32]. 

The polynomial (15) is depicted as dashed line in 
Fig. 4. The solid lines correspond to four different Dlog- 
Pade approximants [55] of this quantity and constitute 
a reliable extrapolation. The plain scries result can be 
trusted up to a; « 0.7. Since each rung can be in four 
different states {s,t^,t^,t~^) the Hilbert space has di- 
mension 4^^ = 2'^*'. Thus the computer calculations used 
about 1 GByte. They took about 20h. The polynomial 
will be made available on our home pages [50]. 



B. One-Triplon Dispersion: Hi 

We define \i) to denote the eigen state of H± with one 
triplon on rung i and singlets on all other rungs. The 
magnetic quantum number m of the triplon at rung i is 
of no importance in the following considerations, since 
HcH conserves m and the total spin S (cf. Tab. I). Thus 
it is not denoted explicitly. 
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FIG. 3. Symbols as in Fig. 2; here: rung-singlets (-triplons) 
are denoted by open (filled) circles. In second order only 
T-2T2 and in third order only T-2T0T2 contribute to eo- The 
action of these operators on the ground state (always to the 
left) is shown step by step. The upper part shows that in 
second order the ground state energy per site becomes inde- 
pendent of the cluster-size, if the cluster contains more than 
two rungs. The lower part shows that wrap-arounds are pos- 
sible, if the cluster is too small: For the third order contribu- 
tion to £() a cluster containing three rungs is undersized. The 
To operator in the middle of the process can break up one 
triplon-pair by moving one of the triplons away by one rung. 
On the three-rung cluster it joins back the remaining triplon 
from the other side (wrap-around) , which would not be possi- 
ble in the thermodynamic limit. This process can be avoided 
by adding one extra rung as buffer, as depicted in the bot- 
tom process. Note that r_2 destroys two triplons only if they 
are neighbours. Thus states with diagonally arranged triplons 
do not contribute. Remembering that one always ends with 
the ground state |0) if one starts from |0) (Hcs conserves the 
number of triplons!), it is clear that this argument works for 
all possible processes T(m) in Hes. 
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FIG. 4. Ground state energy per spin as function of x. 
The plain series result (15) is depicted as dashed line. Four 
different Dlog-Pade approximants ([7,6], [8,5], [5,8] and [6,7]) 
are shown as solid lines. 



on neighbouring sites (triplon-pair) if both of them are 
occupied by singlets. Suppose that T2 is immediately fol- 
lowed by the T_2 operator. Then there can be a hopping 
of the initial triplon by two rungs, if the triplon-pair was 
created in the immediate vicinity of the triplon at site 
i to produce a three-triplon state. The situation is de- 
picted in Fig. 5a. This is a second order process. It 
moves the triplon by two rungs. We could go on like this 
(. . . T^2T2T-2T2) or we could start to build up a linked 
chain by iterative application of T2 operators, say, to the 
right of the triplon at site i and then destroy the chain 
from the left (e.g. T_27^-27272)- All these processes lead 
to a maximum motion of the initial triplon by I rungs in 
Z*^ order. 
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Since Hcs{x) conserves the number of triplons the ac- 
tion of Hcfi{x) on the state \i) is a hopping of the triplon. 
We define the hopping coefficients 



4'^.(x) = mMx)\3) 



(16) 



The superscript cl indicates that the hopping coefficient 
might depend on the cluster on which it was calculated. 

The hopping coefficients ti-j of the irreducible one- 
particle operator Hi read (see Eqs.9 in Ref. [12]) 



ti,j = {i\Hi\j) = {i\H,s - Ho\j} = 4',. - E'o% 



(17) 



Since Hi is a cluster additive, i.e., an extensive, opera- 
tor, the coefficients U-j can be calculated for the infinite 
system on finite clusters up to some finite order. This is 
the reason why we dropped the superscript cl from ti-j. 
The cluster ground state energy Eq^ must be calculated 
on the same cluster as the "raw" hopping coefficients a^}j . 

For each order of the coefficient ti-^j there exists a min- 
imum cluster which must contain the two rungs i and j. 
To classify the size of the minimum cluster we study how 
far the triplon motion extends in a given order I. Only 
processes, which take place on linked clusters of active 
bonds (see previous section), contribute to the extensive 
thermodynamic hopping coefficients tj-j. The minimum 
cluster must be a linked cluster, which contains the rungs 
i and j. 

The action of a single Tq operator (first order process) 
on \i) is to shift the triplon by one as can be readily 
seen from Tab. I. Somewhat more intricate is the case of 
the operator T2 acting on In any operator-product 
T(m) an operator T2 is always accompanied by a destruc- 
tion operator T_2. The operator T2 creates two triplons 
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FIG. 5. Processes of H^g that lead to a motion of the ini- 
tial triplon on rung i. Active bonds are depicted by thick 
lines. All processes that contribute to thermodynamic exten- 
sive hopping coefficients take place on linked clusters of active 
bonds. Part a) shows a second order process moving the ini- 
tial triplon by two rungs. Paxt b) is a process of order n + 2 
moving the triplon by n -|- 2 rungs. 



The creation of a triplon-pair not connected to the ini- 
tial triplon on site i does not lead to any motion of the 
latter unless there is a sufficient number of Tq operators 
moving the triplon at site i towards the isolated triplon- 
pair until they form a state with three adjacent triplons 
as depicted in Fig. 5b. This also leads to a maximum 
motion of the initial triplon by I rungs in Z'^ order. 

All possible combinations of the T2, Tq and T_2 oper- 
ators that can appear in a T(m)-product of f^eff can now 
be viewed as a product of the processes discussed. So we 
conclude 



maximum motion of one triplon under the action 



of Hcs in order = I sites 



(18) 



Therefore, the minimum cluster to calculate the hopping 
coefficient ti-j in order I in the thermodynamic limit must 

contain the two rungs i and j . which must not be further 
apart than / rungs. Additionally the minimum cluster 
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must contain all I bonds that can be activated in all pro- 
cesses moving the triplon from rung i to rung j. Fig. 6 
illustrates the situation for all coefficients that can be 
calculated in fourth order. 
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FIG. 6. AH possible hopping coefficients that can be calcu- 
lated in 4"" order. Again, active bonds are depicted by thick 
lines. All processes that have to be considered take place on 
linked clusters. The initial (final) triplon positions are de- 
picted by a filled circle (cross). They are contained in the 
minimum cluster (cl), which is defined by all active bonds 
for each coefficient. The exclamation mark next to the af^ 
cluster is to remind us that we have to subtract the cluster 
energy Sg' to get the cluster independent hopping coefficient 
U;i=to, c.f. Eq. (17) 



Because of translational invariance of the original un- 
derlying model we can choose a suitable origin on each 
minimum cluster and it suffices to use a single label for 



the hopping coefficients, i.e., ti 



t, 



-1. Addition- 



ally, one has td — t^d due to inversion symmetry. Note 
that these relations follow only for the thermodynamic 
hopping coefficients. In general, the cluster specific coef- 
ficients (t, aY^ have lower symmetries. Following the ar- 
gument above we calculate all thermodynamic hopping 
coefficients in 14*'^ order (tp, ^i, • ■ • ^14) on an open cluster 
of 15 rungs. 

From the thermodynamic cluster-independent hopping 
coefficients we construct the one-triplon energies. We 
define the Fourier-transformed states 



\k) 



1 



(19) 



where j counts the rungs and N is the total number of 
rungs. Calculating the action of Hi on these states yields 



In 



E ^"^^^;4E^-"^'b-) 



(20a) 
(20b) 
(20c) 



|fe) 



Making use of the inversion symmetry td — t^d yields 
the real one-triplon dispersion 



(21) 



uj{k]x) = {k\Hi{x)\k) = to + 2 E tdCOs{dk) , 



where the maximum order Zmax is 14 in our case. 

Again, the hopping coefficients a.i-j and the corre- 
sponding cluster ground state energy Eq and therefore 
the cluster-independent coefficients td are calculated by 
implementing the action of Hcs on the states \i) on a 
computer (details in Ref [7]). For fixed one-triplon mo- 
mentum k the one-triplon dispersion is a 14"^ order poly- 
nomial in X with real coefficients. Thereby we retrieve 
and extend the numerical 13* order result in Ref. [32]. 

In Fig. 7 the one-triplon dispersion is displayed for 
five different x-values. For x —0.2, 0.4 and 0.6 the grey 



2.0 



1.6 



2^ 



1.2 



0.8 



0.4 



' 1 ' 1 ' 1 — 


' 1 ' 






























■ . ■ . ■ 





0.0 



0.2 



0.4 0.6 



0.8 



1.0 



FIG. 7. One-triplon dispersion for various a;-values as indi- 
cated. The grey dashed lines correspond to plain series results 
{x =0.2, 0.4, 0.6) or to optimized series results {x =0.8, 1.0). 
The solid curves depict our most reliable results obtained by 
the novel extrapolation scheme explained in Ref. [57]. 



dashed curves represent the results obtained by using the 
plain series results for the hopping coefficients td- The 
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grey dashed curves for x =0.8 and 1.0 result from using 
the optimized hopping coefficients according to the opti- 
mized perturbation theory explained in Sect. V (param- 
eter choice: Qopt = 2.9a;). These results are compared to 
the most reliable extrapolations depicted as solid lines. 
The latter are obtained by using the novel extrapolation 
scheme based on a re-expansion of the original series re- 
sults (polynomials in x) in terms of a suitable internal 
variable p{x) [57] . 

C. Two-Triplon Interaction: H2 

We define the states denoting the eigen state of 

H±^ with triplon 1 on rung i, triplon 2 on rung j and 
singlets on all other rungs. Two triplons together can 
form an 5 = singlet, an S* = 1 triplet or an 5 = 2 
quintuplet bound state. Tab. II summarizes these nine 
states sorted by their total spin S and magnetic quantum 
number m. 

By construction iJgff conserves the total spin S and the 
magnetic quantum number m. Therefore it is convenient 

to work in the basis given in Tab. II. This table defines 
the states by the linear combinations in the third 

column. 



s 


m 




2 
2 
2 
2 
2 


2 
1 

-1 

-2 


l/V2(|t«,t^) + |i\t«)) 
l/^{\t-\t^)+2\t\t') + \t\t~^)) 
l/^[\t-\t°) + \t\t-')) 
\t-\t-^) 


1 
1 

1 


1 


-1 


l/^2{\t\t^)-\t\t')) 

l/^/2(|^^^-l)-|^-^^l» 

l/^(|t°,i-i)-lt-i,iO) 








i/Vi{\t'',t'')^\t\t-^)-\t-\t^)) 



TABLE II. The nine states two triplons can form, com- 
bined to states with given quantum numbers S and m. 



Again, due to triplon conservation the action of ifeff 
on the state is to shift the triplons to rung i' and 
rung j' conserving also S and m. Nothing else is possible. 
In analogy to Eq. (16) of the preceding section we define 
the interaction coefficients 

o^m^= '{^,3\HM\k^f ■ (22) 

The coefficients depend on the total spin S but not on 
the magnetic quantum number m. Hence the m-index is 
dropped here and in the following. 

The exchange parity is determined by the total spin S 

\i,jf = {-lf\j,if. (23) 

This means that we can restrict the description to those 
states |i, j) for which i < j. 



Making use of the above the irreducible two-triplon 
interaction coefficients tfj.f^i follow from 

tfr,ki = ''{i,j\H2\k,lf = %j\H,fi -H,- Ho\k,lf 

S,c\ zT'cl r r 
=%;kl - ^0 Oi,kOj-l 

(24) 

analogous to Eq. (9) in Ref. [12]. Again, Eq^ and the 
one-triplon hopping coefficients tj.j must be calculated 
on the same cluster as the "raw" two-triplon coefficients 
^fj-ki- The cluster hopping coefficients t^li are needed 
only in the intermediate steps of the calculation of the 
irreducible interaction coefficients. 

There will be no tu-ki or tij-kk since it is not possible 
to have two triplons on one rung at the same time. This 
constraint can be viewed as a hardcore repulsion interac- 
tion. 

The construction of the minimum cluster needed to cal- 
culate the Uj-ki in the thermodynamic limit follows the 
same line of argumentation as in the one-particle sec- 
tion. Generally, the cluster must be large enough to en- 
compass all possible processes in order I. The minimum 
cluster has to include all linked bonds that can bo acti- 
vated in any possible interaction process of length I which 
leads to state if one starts with state Obvi- 

ously the rungs i, j, i' and j' must be contained in the 
minimum cluster and they must be connected by active 
bonds. Fig. 8 shows some interaction coefficients with 
fixed initial configuration (adjacent triplons) and their 
associated minimum clusters in 4**^ order. All interac- 
tion coefficients of order I can be calculated on a cluster 
containing I + 1 rungs. 

For particular systems there may be symmetries, e.g., 
spin rotation invariance, or other particularities, e.g., 
nearest-neighbour coupling only, which prevent certain 
processes from generating non-vanishing coefficients. In- 
deed, we found that the spin ladder with nearest- 
neighbour coupling in order I induces only interaction 
coefficients of a range which can be determined by dis- 
tributing I hops among the two triplons. For instance, 
in 4**^ order the triplon at rung i may hop to rung i — 1 
and the triplon at rung j hops to rung j + 3. Another 
possible process might be that the triplon at rung i stays 
at this rung while the triplon at rung j hops to rung j -f- 4 
or maybe only to rung j + 3 and so on. This particular- 
ity implies for instance that the process shown in Fig. 8c 
vanishes for the spin ladder system. 

Due to the translational invariance of the ladder sys- 
tem the momentum A; is a good quantum number in the 
one-triplon sector and the diagonal matrix elements of 
the Fourier transformed states jfc) are the cigcn energies 
u>{k). With two triplons present only the total momen- 
tum K is & good quantum number. The relative mo- 
mentum q is not conserved and generally leads to the 
formation of a two-triplon continuum. 

To make use of the conserved total momentum we turn 
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FIG. 8. Some interaction coefRcients that can be calculated 
in 4**^ order. Initial (final) triplon pairs are denoted by full cir- 
cles (crosses). The thick solid lines are active bonds between 
rungs. Exclamation marks indicate that one has to subtract 
one- or zero-triplon terms according to Eq. (24) to obtain the 
extensive thermodynamic interaction coefBcients. 



to a new basis. As a first step we use center-of-mass 
coordinates, i.e., \r,r + d)^ = {—l)^\r + d, r)^ , 

with r — i and d — j — i. The restriction i < j (see text 
below Eq. (23)) translates to rf > 0. We choose a suitable 
origin, say k — 0, and rename 



= {r,r + d'\H2\0,d) = {i,j\H2\k,l) = Uj-m 



with d — I, r = i and r + d' = j. From Eq. (24) we then 
obtain 

id;r,d' = '^d-f,d' ^ £'o''^0,r<^(i,r+d' — tf.QSd,r+d' — t'd;r+d'^0,r 
- 4,r+d'SdA-^f - tiMr+d'{-lf (26) 

in the new basis. This is equivalent to the equations 
emerging from considering the special cases 



^d;0,d' — "'d;0,d' ~ ''d'-d ~ "d^d'y^Q ' ^ I 



(27a) 



td:d-d',d' = o-d^d-d'.d' ~ ^d-d' - ^d,d'{tt + E''^) (27b) 



t 



d-.~d' ,d' 



S,c\ _ ,cl 
"'d-.-d'.d' ''-d-d' 



(-1)' 



t 



d-d.d' 



S,c\ _ 
'^d;d.d' ''d+d 



ttd'i-^r 



(27c) 
(27d) 



Otherwise the interaction coefficients t^.^. ^, and a^.^ ^, are 
identical. 

As a second step the states \r,r + d)^ are Fourier 
transformed with respect to the center-of-mass variable 
(r + d/2) 



(-1) 

r^r-\-d 



1 



N ^ 
1 

'N 

= (-l)^|X,-d)' 



\r + d,rY 

iK(r-d/2) I 



(28) 



where K is the conserved total momentum in the Bril- 
louin zone and N is the number of rungs. For fixed K 
and S the relative distance d > between two triplons 
is the only remaining quantum number one has to keep 
track of. 

To obtain the complete two-triplon excitation energies 
we have to calculate the action of 



Hr 



Ho — Hi + H2 



(29) 



on the two-triplon states \K,d). The two addends on the 
right hand site are considered separately in the following. 

The operator Hi can move one of the two triplons at 
maximum. A short calculation yields 



Hi\K,d)' 



.iKir+d/2) 



tni\r + n,r + d)'=^ + \r,r + d-7i)^) 



/2 I ^~iKn/2 



1 ye»i^('-+(d-n)/2)| 



K,d—n) 



(25) ^ 



tmax 

2 £ t„cos(if|) [sgn(d-7i)]^|if,|d-n|)« . (30) 



Here we used the previously calculated matrix-elements 
tn — t-n (inversion symmetry), which have been cal- 
culated to Imax = 14 (cf. preceding subsection). Since 
we restricted d > the sgn-function enters the result 
by Eq. (28). For fixed Hi now appears as a semi- 
infinite band matrix in the remaining quantum number 
d. Independent of the size of the initial distance d > be- 
tween the two triplons, Hi will produce states where the 
distances between the triplons are incremented or decre- 
mented by 14 rungs (^max — 14) at maximum. If the 
initial distance d is larger than 14, Hi continues to pro- 
duce the same matrix elements on and on for all d > 14, 
i.e., the matrix representing Hi in the chosen basis for 
fixed K is semi-infinite with a repeated pattern in the 
tail. The head of -ff 1 , i.e., the 14 x 14 block between states 
with d < 14, contains matrix elements with a somewhat 
more complicated structure. Here the matrix element 
between the starting distance d and the final distance d' 
is a sum of the direct process d ^ d' , where one of the 
triplons has hopped n rungs to the right (n > 0) or to 
the left (n < 0) with d — n = d' > 0, and the indirect 
process with d — n — ~d' < 0. The situation is sketched 
in Fig. 9. The matrix Hi comprises the full thermody- 
namic one-triplon dynamics in the two-triplon sector for 
the given order I max = 14. 

The situation is more complex for H2. In a first step 
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we find 



H2\K,df = 



max { n + d' , d — n } 
<!max 



(31) 



max{ n+d',d — n} 
^^max 



with the two integers n G Z and d' G N as summation 
indices. The positive distances d and d' must be smaller 
or equal to Imax, since a maximum of Imux linked bonds 
can bo produced in this order and all four triplons sites 
(the two initial sites and the two final sites) must be 
contained in the resulting linked cluster. The td-n,d' are 
the matrix elements of H2 defined in Eq. (26). The last 
equality follows from substituting the summation-index 
r ^ r + n. 

To simplify the expression further inversion symmetry 
is used. We have 



td-y,d' = {r",r" +d'\H2\r,r + d) , 



(32) 



with r' = r" — r. The thermodynamic interaction coef- 
ficient td;r',d' is associated with a fixed constellation of 
initial and final triplon pairs. We define a configuration 
CON by the set of four positions given by these two pairs 
CON = {r,r + d,r",r" + d'}. Let s denote the middle 
of this configuration s = (max(CON) — min(C0N))/2. 
Reflecting a configuration about s and interchanging the 
triplon positions in both initial and final triplon pairs 
gives 



td;r',d' = {r",r" + d'\H2\r,r + d) 

={2s - r" - d', 2s - r"\H2\2s - r ■ 

= td-d-d'-r' ,d' ■ 



d,2s- 



(33) 



Possible minus signs cancel since they appear twice. We 
can now split the sum over n in Eq. (31) in three parts, 
n> {d- d')/2, n<{d- d')/2 and n={d- d')/2. The 

second sum is indexed back to n > (rf — d') /2 by making 



use of J2„ 



Y.n>i "2j-n where j := {d - d')/2 



H2\K,df=J2 



x:{n+d',d-n}<imax 
n>(d-d0/2e2 



+td;d-d'-n,d'e'''^''-^''-''^/'^\K,d'f 

td;(d-d')/2,d'\^,d')^ 



nax{n+<i',<i-n}<!„ 
n=(<i-<i')/2eZ 



= 2 J2 td;n,d' cos[K{n -{d- d')/2)]\K, d') 



<i\S 



nax{n+<i',<i-»i}<!„ 
n>(<i-<i')/2£Z 



+ X/ ^d;{d-d')/2,d'\K,d')^ . 



(34) 



nax{n+d',(i — n}<iniax 
n = ((J-(i')/2eZ 



In contrast to Hi the matrix representing H2 is of fi- 
nite dimension due to the finite range of the contributing 
processes (finite maximum order) expressed by the re- 
strictions of the sums appearing in Eq. (34) . In our case 
the td;r,d' have been calculated up to Zmax = 13 giving 
rise to a 13 X 13 matrix in the distance d for fixed K. 
Fig. 9 sketches the situation. 

Finally, the sum of the two matrixes H\ and H2 with 
respect to basis (28) comprises the complete two-triplon 
dynamics. 

The above approach is well justified. At large dis- 
tances the two-triplon dynamics is governed by indepen- 
dent one-triplon hopping. At smaller distances an ad- 
ditional two-particle interaction occurs given by td-r,d' 
connecting the state \r,r + d') with state \Q,d). The sum 
Hi + H2 gives the combined effect of one-triplon hopping 
and two-triplon interaction. Fig. 10 shows that the in- 
teraction coefficients {K, d.\H2\K, d) for the ladder system 
indeed drop off rapidly for larger distances. 



head 








\ / 



FIG. 9. The left part of the figure schematically shows 
the matrix representation of H\ and H2 in the two-triplon 
{\K, d)} basis (28). The matrix H\ has elements in the whole 
grey area, while H2 has elements in the dark grey area only. 
We calculated the elements of H2 up to order 13 and those 
of H\ up to order 14, so that H2 is a finite 13 x 13 matrix 
and H\ a semi-infinite band matrix, whose width is 28, see 
Eqs. (30) and (34) for further information. The sum of 
and H2 represents Hcs in the two-triplon sector to the given 
orders. The right part shows the initial vector |Init) = Oeff |0) 
as calculated in Sect. Ill for the two-triplon sector. Since we 
calculated OcB up to order 10, |Init) is a vector of dimension 
10 in the {\K, d)} basis. The Grccu function Q (Eq. (3)) is cal- 
culated by tridiagonalizatiou, more information in Sect. IV. 
For K and x fixed, the elements of the matrix and the vector 
reduce to real numbers. 



Taking the perturbation expansion up to order Zmax 
allows to calculate the irreducible two-particle interac- 
tion up to a distance I between the two particles cor- 
rectly within order Zmax- No processes involving larger 
distances appear. But the part of the two-particle sector 
that can be described by one-particle dynamics alone is 
taken into account for all distances between the two par- 
ticles and describes hopping processes of range < ^max 
correctly within order ^max- 

Pit K = 11, the smallest eigen value of Hi+ H2, where 
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FIG. 10. Expectation values of some diagonal elements of 

H2 (pure two-triplon interaction) for the ladder system in 
the {\K,d)} basis for various if-valucs and x set to 0.6 as 
function of the remaining quantum number d. Clearly, the 
irreducible two-triplon interaction coefficients drop off rapidly 
with increasing relative distance d between the two triplons. 
Non-diagonal elements behave in a similar fashion. At larger 
distances two triplons are asymptotically free. 



Hi is the upper left 13 x 13 sub-matrix of Hi, can be 
extracted as a 13**^ order polynomial in x and identified 
as the (lowest) bound state. This is possible because at 
K = n, the relative motion of the two triplons is of order 
while the interaction enters in order x. Hence the 
interaction dominates over the kinetics for a; ^ so that 
a local bound pair is the simple eigen state for vanishing 
X. Our results extend the results by Zheng ct al. [11] for 
S = 1 from 12*^^ to 13*^^ order, and for 5 = from 7*'^ to 
13*^ order. The polynomials will be made available on 
our web pages [50]. 

III. TRANSFORMATION OF THE 
OBSERVABLES 

A. General Aspects 

The continuous unitary transformation of observables 
has been explained in detail in Rcf. [12]. Here we briefly 
review the most important general aspects before we de- 
scribe the procedure for the spin ladder in detail. 

Using the same transformation as for the Hamiltonian 
we derive a series expansion (similar to Eq. (11)) for the 
effective observable Oca onto which a given initial observ- 
able O is mapped by the perturbative CUT procedure 

OMx) = J2^''J2 E C{m;i)0{rn;i) , 

k=0 i=l |m|=fe (35) 



0{m; i) := • • • T„,_,C>T„, • • • . (36) 

The operators Tj are the same as in the Hamilton trans- 
formation. The coefficients C(m; i) can again be calcu- 
lated on a computer. They are fractions of integers [12]. 

The effective observables are described by weighted vir- 
tual excitation processes T{rn) interrupted by processes 
induced by the observable as given in Eq. (36). Some- 
times it is convenient to seek for a decomposition of O 
with respect to its action on the number of particles 

AT 

0= J2 ^n, (37) 
n=-N 

where creates n particles or destroys them if n < 0. 

An important point is that Ocff is not an energy quanta 
conserving quantity, i.e., it does not conserve the number 
of triplons in the spin ladder system. This is formally ex- 
pressed by the fact that the sum over m is not restricted 
to M(m) = 0, so that OeS can add or subtract an arbi- 
trary number of particles. 

The effective operators OeS can be decomposed in a 
sum of cluster-additive operators Op,n, for which the 
linked cluster theorem can be used 

00 

n=0 p>—n 

Here p indicates how many particles are created {p > 0) 
or destroyed {p < 0) by Op,n- The subindex n > in- 
dicates the minimum number of particles that must be 
present for Op,„ to have a non zero action. The action of 
the operator Op,n on a state containing less than n par- 
ticles is zero. Further definitions and details concerning 
the operators Op^n can be found in Ref. [12]. 

Focusing on T = experiments in the following, we 
treat only the operators C'p>o,o- Their interpretation is 
particularly simple. The effective observable Ocff act- 
ing on the T = state, i.e., the ground state or excita- 
tion vacuum, respectively, decomposes in a sum of the 
operators Op>o.O! each injecting p = 0,1,2... triplons 
into the system. In Ref. [12] we showed that the C'p>o,o 
can be directly calculated from the action of Oes on |0) 
on minimum finite clusters. No extra terms have to be 
subtracted to obtain thermodynamic results. The cal- 
culations can again be implemented on a computer in 
analogy to what was done for the Hamiltonian. 

To be more specific let O be a locally acting observable 
injecting triplons at a specific site r of the ladder. Then 
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the effective observable reads 

Oeff(r)|0)=^Op,o(r)|0) 

p>0 



= c|0)+ Cn\r + n) + 

+ X] Cn,n'\r + n,r + n') + ■ 



|n| + |«'|<!„ 



(39) 



The restriction |n| + < l^^x for the third sum reflects 
the fact that the two triplons, after being injected, can- 
not undergo more rung-to-rung hops in total than the 
maximum order /max- Therefore, the maximum distance 
p = \n — n'\ occurring is I in order I for the spin ladder 
system. 

Once the coefficients c are calculated the spectral 
weights In are accessible, which are contained in the dif- 
ferent particle-sectors characterized by the number N of 
particles injected 

/jv = (0|0-jv,o(r)Oiv,o(r)|0) 

= J2 Kr- + ni,... ,r-f-njv|Ojv,o(r)|0)|2 



ni ,. . . ,njv 



(40) 



ni,... ,nAr 



If the total weight /tot of the operator is also known, for 
instance via the sum rule /tot = {0\O^\0) - {0\O\0y , the 
relative weights of the individual particle sectors In / /tot 
can be calculated. They serve as an important criterion 
to judge the applicability of our approach. If most of 
the weight can be found in sectors of low quasi-particle 
number and sectors of higher particle number can be 
safely neglected the approach will work fine. The cho- 
sen particles constitute a suitable basis to describe the 
system. This argument has been used by Schmidt and 
Uhrig [18] to show, that the triplon is a well suited par- 
ticle to describe the one- dimensional spin chain. There 
basically all the spectral weight is captured by one and 
two triplons. 

So far local observables 0{r) were considered. A real 
experiment, however, couples to the system in a global 
fashion. Due to translational invariance the injected par- 
ticles (here triplons) have a total momentum K. Thus 
we define the global observables in momentum space rep- 
resentations 

Oeff(if)|0)=^Op,„(/^)|0) 



p>0 



N 



= E-^E^^''^^P.o(r)|0), 

p>0V^r=l (41) 

where A'' is the number of rungs in the system. Let us 
investigate the one- and two-triplon sectors separately. In 



the one-triplon sector we have (here K is the one-triplon 
momentum k) 

Oi,o(fc)|0) = -^5^e''='-^c„|r + n) 

r n 
n ^ r 

= ^c„e-^'="|fc) . (42) 

n 

We used the same definition for |fc) as introduced in 
Sect. II. Due to inversion symmetry c„ = c_„ holds. 
Thus Eq. (42) simplifies to 



(fc|C'i,o(A;)|0) ^Ak = co-H2^c„cos(/cn) . 

n 

Somewhat more complex is the two-triplon sector 
1 ^ 

O2,o(^)|0) = ^y2e''"-y2cn,n'\r + n,r + n') 
ViV — ; 

r— 1 n.n' 



(43) 



■iK{n+d/2) _ 



^E' 



iK{r+d/2) 



n,d 

= YAK,d\K,d) , 



r,r + d) 



(44) 



where we defined the relative distance d = n' — n between 
the two injected triplons. The definition of \K, d) is taken 
from Sect. II. Again, inversion symmetry, here c„_„' = 
(— l)'^c_„,_„', can be used to obtain real results for the 
coefficients AK,d- The variable S G {0,1,2}, which is 
a good quantum number, denotes the total spin of the 
injected triplon pair. 

The action of OeS from the ground state into the two- 
triplon space produces the states \K, d) with < d< ^max 
in order /max- Thus, for fixed K, the action of OeS may be 
visualized as a vector in the remaining quantum number 
d of which the first imax entries are the AK,d of Eq. (44). 
All other entries are zero. This vector, labeled initial 
vector jlnit) for reasons given in Sect. IV, is depicted in 
Fig. 9 together with the matrix representing H^g for fixed 
K in the two-triplon sector. 



B. Observables in the Spin Ladder 

We now turn to the evaluation of the observables of 
interest in the ladder system. We calculated the C{rn;i) 
in Eq. (35) up to and including order /max = 10 for the 
problem under study. The four local operators considered 
are 
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0\r) = Si,,S2,. = 



1IV 



(45a) 
(45b) 

(45c) 
(45d) 



where the decompositions are either given in Table I for 
the r or in Table II for the T^, with n G {I, II, III, IV}. 
The index ? = 1, 2 in Eq. (45b) denotes the leg on which 
the observable operates. 





\s) 
\V) 




-3|s) 


ATI' 


\sX) 
\t\s) 




\t',t^)-\t\t') 
-\t\t^^) + \t\t^) 


\t\s) 
\s,t-^) 




\t\t-') + \t-\t^) 
-\t\t-^) + \t-\t^) 
\t\t-')-\t-\t^) 


\t-\s) 













-'0 
















TABLE III. Action of the local operators T^^ appearing in 
Eqs. (45). The notation is the same as in Table I. 

We start with a simple symmetry property. Let V 
denote the operator of reflection about the center-line of 
the ladder as depicted in Fig. 11. If |n) denotes a state 



FIG. 11. The operator V reflects about the depicted axis. 
A single rung-singlet (-triplon) has odd (even) parity with 
respect to V. The action of V on the rung-singlet ground 
state is defined to be of even parity PjO) = |0). If in |0) one 
singlet is substituted by a triplon we get the state |1) and 
V\l) = Generally, one has V\n) = (-l)"|n). 



with n rungs excited to triplons while all other rungs are 
in the singlet state we find V\n) = (— l)"|n), see caption 
of Fig. 11. The state |n) might be a linear combination of 
many n-triplon states so no generality is lost in writing 



(46) 



n>0 



The parity of the ladder observables introduced in 
Eqs. (45) with respect to V is clear from their definition: 



O^^^ is odd while O' and O^^ are even with respect to V, 
just as the symmetriezed observable O^^ = {O^ + 0^)/2. 
These parities are conserved under the CUT so that V 
applied on both sides of Eq. (46) yields 



0„|O\-/ E„|2«), Ocffcvcn 
^" lEj2ri + l), Oeffodd ■ 



(47) 



We thus find that an even (odd) parity of OcS implies 
that Ocff can inject an even (odd) number of triplons 
into the system. 

The coefficients c in Eq. (39) have been calculated for 
the one- and two-triplon case on a computer in a similar 
fashion as the coefficients t for the effective Hamiltonian. 
The implementation of Ceff acting on the ground state 
|0) follows the same line as described in detail for Hes in 
Ref. [7] . The minimum clusters necessary for some fixed 
order arise from the same considerations as in Sect. II. 
Again, the coefficients c are rational numbers which we 
computed up to imax = 10*^ order for the observables in 
Eqs. (45). 

First physically interesting quantities are the spectral 
weights of the observables. As illustrated in Eq. (40) the 
spectral weight contained in the A''-triplon sector, 1^, 
is readily given by the coefficients c. Under certain cir- 
cumstances the total weight Jtot = Iq + h + I2 + ■ ■ ■ of 
an observable might be accessible from sum rules. Let us 
consider the S = operator Ojg- in Eq. (45b) as an exam- 
ple. Here the total weight Itot {x) can be obtained from 
the ground state energy per spin eo(a;) given in Eq. (15). 
Since 20j|f.(x) = d/dxH^six) (cf. Eqs. (1) and (6)) the 
sum rule can be expressed in terms of the effective Hamil- 
tonian, giving rise to 



/tot =}_^In = {0\0' 

N=0 

_ 3 y y2 
~ 16 ~ 2" ~ ^ ' 



{Q\0\Qy 



(48) 



with Y := deo/dx. If both, /tot and some of the In are 
known, we can calculate the corresponding relative spec- 
tral weights In /hot as functions of x. Fig. 12 shows the 
resulting relative weights for the observable Oj^ for the 
first four triplon sectors. Since one cannot form an = 
object from a single rung-triplon there is no Ii for this 
observable. The contribution of I^/Itot is of order 10~^ 
leading to no visible changes in Fig. 12. Contributions of 
higher triplon channels are expected to be even smaller. 
All relative weights add up to unity. As can be seen in 
Fig. 12 the first four relative weights fulfill this require- 
ment with great precision. For x = the singlet made 
from two isolated-rung triplons contains the full weight of 
the considered operator. As x increases the triplons start 
to polarize their environment, the two-triplon weight de- 
creases and multi-triplon states gain weight. A similar 
figure for the S = 1 operator Sf ^ can be found in Ref. [13] 
(Fig. 2). 
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FIG. 12. Relative weights for the S=0 operator S'i,iSi,i+i 
(Eq. (45b)). The /n are calculated according to Eq. (40) up to 
and including order 10, 8 and 7 in a; for the C{„} for =2, 3 
and 4, respectively. The total intensity Jtot has been extracted 
from the 14"^ order result for the ground state energy per spin 
according to Eq. (48). This figure represents an improvement 
of a similar figures in Ref. [13]. The expansion for I3 has been 
extended by one order and the extrapolation of I2 has been 
improved. 



tion of the relative distance d for fixed iiT- values at a; = 1. 
The depicted values are obtained by using standard Pade 
techniques for the AK,d as polynomials in x for fixed mo- 
mentum K. All calculated one-triplon (Ak) and two- 
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FIG. 13. The two-triplon momentum dependent coeffi- 
cients Aji^d of the observables O^^ and O^^ for all calculated 
distances d at a; = 1 and momenta as indicated. The coeffi- 
cients rapidly drop to zero with increasing distance. It is not 
necessary to go to larger orders, i.e., distances. 



From the depicted result we conclude that the triplon is 
an excellent choice for quasi-particle in the ladder system. 
For X not too large most of the spectral weight is captured 
by a few triplons. Calculations containing only a few 
triplons suffice to explain most of the physics for x ^ 1.5. 

Eqs. (43) and (44) show how the momentum depend 
coefficients Ak (one-triplon) and Af^d (two-triplons) can 
be calculated from the corresponding c-coefficients. For 
the two-triplon sector we provide some examples to 3'"^ 
order in x. The shown coefficients Ax^d belong to the 
symmetrized observable O^^ — {O]^ + 0}^)/2 or to the 
observable O^^ , respectively. 



A^^ - 



A 



K,2 



4 IV 



1 1 

X 

4 8 



1 /25 17 3 
8 U + T^'^^^^^'V" 
icosQx).+ lcosQx) x^ 

1 f 37 (1 \ 13 /3 

— cos - A cos - K 

16 V 16 V2 / 16 V2 



1 . /I 

— X sin — K 

2 V 2 



11 3 • n 

— x"^ sm -K 
64 12 



1 2 A 
- X sin - 
4 \ 2 



(49) 
(50) 



triplon {AK,d) coefficients will be made available on our 
home pages [50]. 



IV. EVALUATING THE GREEN FUNCTION 

We are now in the position to calculate the zero tem- 
perature one- and two-triplon spectral densities associ- 
ated with the ladder observables introduced in the last 
section. To this end we start by analyzing the energy and 
momentum resolved retarded zero temperature Green 
function 



■00 



0\K) 



uo - {H{K) - Eo) 



iO- 



-0{K) 



(51) 



from which the spectral density S{K, w) follows by taking 
the negative imaginary part. As explained in Sect. LB, all 
operators can be replaced by their effective counterparts 
after the transformation and the ground state by the 
triplon vacuum 1 0) . 



Fig. 13 displays the coefficients of O^^ and O^^ as func- 
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A. One-Triplon Green Function 

In the one-triplon case the calculation is particularly 
simple. Using Dirac's identity 



= V- 



=F iT^5{x — xq) , 



X — Xq± iO+ X — Xq 



where V denotes Cauchy's principal value, we find 



(52) 



Ol,{k)S{c, - H,)0,^,{k) 

\Akf{k\S{w-H,)\k) 
\Ak\''S{uj-Lu(k)) . 



(53) 



The one-triplon dispersion w(fc) and the observable coef- 
ficient Ak are readily given by Eqs. (21) and (43), respec- 
tively. At each point {k,uj{k)) the corresponding weight 
is given by the square of the modulus of Ak which is a 
polynomial in x. The result is thus obtained by assign- 
ing a (5-function with corresponding weight to each point 
{k,u{k)). 



B. Two-Triplon Green Function 

For the two-triplon case we choose to evaluate the effec- 
tive Green function by tridiagonalization. This leads to 
the continued fraction expression ( [58-60] , for overviews 
see Refs. [61,62]) 



{G\0UK)0^JKm 



u) — ao 



u) — ai — 



hi 



(54) 



uj — ■■ ■ 



where we can also write for the expression in 

the numerator on the right hand side. The coefficients 
-^K,d are given by Eq. (44). The coefficients and are 
calculated by repeated application of Hqs—Hq = H1+H2 
on the initial two-triplon momentum state jlnit) = |/o) = 
C2,o(^'')|0). The action of Hi and H2 on these states 
have been calculated previously. The results are given in 
Eqs. (30) and (34), respectively. 

Setting the state |/-i) = the recursion (Lanczos 
tridiagonalization) 

= {Hi + H2)\fn) - anlfn) - |/n-l) , n G N 

(55) 

generates a set of orthogonal states if the coefficients 
are defined according to 



an = 



{fnl iHi + H2)\fn} 
iUfn) 



bl+i = 



(/rt+ll/n-n) 
iUfn) 



In the generated {|/„)}-basis H^ff is a tridiagonal matrix, 
where the a, are the diagonal matrix elements and the 



bi are the elements on the second diagonal. All other 
matrix elements arc zero. 

Fig. 9 illustrates the procedure for the two-triplon sec- 
tor. For fixed K the relative (positive) distance d between 
the two injected triplons is the only remaining quantum 
number. In this basis Hi + H2 is represented as a matrix 
(left side). The matrix elements are polynomials in the 
perturbation parameter x. We have to apply this matrix 
iteratively to the !/„). The components Aka of the ini- 
tial vector 1/0) = |Init) are polynomials in x for fixed K. 
By this procedure a new basis |/„) is generated in which 
the fairly complicated matrix in Fig. 9 is simplified to a 
tridiagonal one. 

The general case of more than two triplons can be 
treated similarly. For n triplons we have to consider 
the conserved total momentum K and n — 1 relative dis- 
tances. Then, for fixed K , |Init) and iJgff are still repre- 
sented by a vector and matrices, respectively. But they 
are more complicated. For three triplons we obviously 
have to apply Hi + H2 + H^ to |Init). For four triplons 
Hii is added and so on. The calculation of if„ with n > 2 
is indicated in Ref. [12]. 

Thus, the full many-particle problem is effectively re- 
duced to a few-particle problem! Further, after fixing the 
parameter value x the coefficients ai and hi are obtained 
by the numerical Lanczos tridiagonalization, which is left 
to the computer. This procedure is realized for fixed x 
and K. For the spin ladder we were able to implement 
a maximum relative distance d of « 10000, allowing to 
repeat the recursion about 650 times giving the first 650 
coefficients and 6f . 

Thus the chosen method to evaluate the effective Green 
function introduces no quantitative finite size effects. 
The problem of calculating the spectral densities for 
given effective Hamiltonians and observables comprises 
the two quantum numbers K and d in the two-triplon 
sector. For each triplon more, there is one more relative 
distance to be considered, see above. Our calculations 
are exact and without finite-size effect to the given order 
for the total momentum K. 

There are two approximations that involve the rela- 
tive distance d. They will be discussed in the following. 
The prevailing approximation is caused by the finiteness 
of the perturbative calculations. The true many-triplon 
interactions are accounted for only if all involved par- 
ticles are within a certain finite distance to each other. 
This approximation is controlled, since one generically 
observes a rather sharp drop of the interaction matrix el- 
ements with increasing distances (see Fig. 10 for the spin 
ladder). Gapped systems with finite correlation lengths 
are well suited to be tackled by our method. Difficulties 
arise if the correlations drop slowly with increasing dis- 
tances. In this case the truncations in real-space might 
not be justified. 

A second minor approximation to the results involv- 
ing the relative distance d is introduced by truncating 
the continued fraction expansion of the Green function. 
However, allowing for distances of up to 10000 lattice 



15 



spacings as in the ladder example guarantees that this 
error is extremely small in comparison to the error intro- 
duced by truncating the perturbative expansion as dis- 
cussed in the preceding paragraph. 

The finiteness of the continued fraction can be partly 
compensated by suitable terminations as shall be ex- 
plained in the following subsection. 



C. Terminator 

The spectral density S{K,uj) at fixed K as obtained 
from the truncated continued fraction (54) of the effec- 
tive Green function has poles at the zeros of the denom- 
inator. Thus, S would be a collection of sharp peaks. 
A slight broadening of S via co ^ to + iS {5 small) in Q 
will smear out all poles to give a continuous function for 
all practical purposes. However, we can achieve perfect 
resolution of S as continuous function by introducing a 
proper termination of the continued fraction exploiting 
the one-dimensionality of the considered model. 

For fixed total momentum K the (upper) lower band 
edges (cub) eib of the two-triplon continuum can be calcu- 
lated from the one-triplon dispersion wi (21). All energies 
of the two-triplon continuum are seized by 

cj2iK, q) = wi {K/2 + g) + {K/2 - q) , 

(56) 

where q G [— tt, tt] denotes the relative momentum. 
Therefore, we can calculate eub and eib from the one- 
triplon dispersion 



and define 



euh{K) = m.ax{uJ2{K,q)) 

Q 

eih{K) = mm{u2{K,q)) . 

Q 



(57) 



For fixed K the upper and lower band edges e^b and 

eib determine the values to which the continued fraction 
coefficients and bi converge for i ^ oo. One finds a^o = 
(cub + eib)/2 and 6oo = (cub - eib)/4 [61]. This serves as 
an independent check for the calculated coefficients. 

If we assume the system under study to be gapped the 
massive elementary excitations show quadratic behaviour 
at the dispersion extrema. Then it is generic that u)i {q) is 
smooth, i.e., two-fold continuously differentiable, and so 
is uj2{K^q)- Since the problem is one-dimensional there 
are square-root singularities in the density of states at the 
edges of the continuum if the two particles are asymp- 
totically free at large distances. In conclusion, a square 
root termination for the continued fraction is appropri- 
ate [61,62]. The listed properties lead to a convergent 
behaviour of and 6? with 



ai 
hi 



a„o + 0(lA^) 

&oo + 0(lA^) 



(58) 

and it is well justified to assume ai and hi to be constant 
beyond a certain (large) fraction depth i. Hence we use 



£) = 46^ - (a; - a. 



(59) 



T = (uj-Goo- \/^) for uj > 

r = (uj -ttoo- iVo^ for eib < w < Cub 

T = -4- - Ooo + \/-d) for CO < eib (60) 



The last calculated 6? in Eq. (54) is multiplied by the 
appropriate terminator r. Taking the imaginary part of 
the resulting expression for the case within the contin- 
uum yields the continuous part of the spectral density 
S in the thermodynamic limit. The result is a continu- 
ous function displaying the full weight of the continuum 
correctly, limited only by the finite order of the series 
expansion. 

In the case of bound states the Green function can be 
written as {K is assumed fixed) 



(010^,002,010) 

w - /(w) 



(61) 



where the function /(w) is real-valued for w < eib. The 
position of possible bound states is given by the zeros of 
g{u)) = u> — f{u>). Let too be such a zero of g. We expand 
g about in w — ojq to first order which is sufficient for 
small deviations from uq 



(O|0log2,olO) 
(u) - wo)(l - d^f{ujvi)) 



(62) 



If Q''-' is the retarded Green function the Dirac-identity 
yields 



«5(w)la-«a 



U .o, , (0|g|,og2,o|0) ,, , 
TT 1 - du;f{u>o) 



(63) 



clarifying that a possible bound state shows up as a 6- 
function. Its spectral weight is given by 



Cund = 5u.(^;^(u.)-^) 



1 - d^f{(Jo) 



(0|0lo02,o|0) '(64) 



which is easy to calculate once the continued fraction is 
known. 

The methods explained in this section have been used 
to derive the spectral densities presented in earlier pub- 
lications; see Refs. [17,13,15,14,19,20,16,21]. Finally we 
address the necessary extrapolations if the perturbation 
parameter x is not small. 



V. OPTIMIZED PERTURBATION THEORY 

The results for the one-triplon dispersion ui{k, x) and 
the matrix elements of Hi and H2 in the two-triplon sec- 
tor are perturbative. They rely on effective operators cal- 
culated as truncated series, i.e., polynomials, in x. The 
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theory is controlled in the sense that it is correct for 
a; ^ 0. But in general we do not have information about 
the radius of eonvergenee. A standard way to extrapo- 
late the polynomials is the use of approximants like Pade 
- or Dlog-Pade - approximants and others [55] . This is a 
feasible task if one is dealing with a few quantities only. 
However, for the matrices Hi and H2 of Sect. II there are 
more than 100 matrix elements, each a truncated series 
in X, to be extrapolated for each K. Clearly, this task 
has to be automatized. The Pade-methods do not allow a 
simple automatization, since the resulting approximants 
are not sufficiently robust. Some of the possible Pade ap- 
proximants to a given polynomial might display spurious 
singularities and there is no way to predict this to hap- 
pen. Pade approximants need to be inspected manually. 

Some progress can be made by a recently developed 
extrapolation procedure introduced in Ref. [57]. This 
technique relies on re-expanding the initially obtain trun- 
cated series expansion in the external perturbation pa- 
rameter X by a suitable internal parameter p{x) . The lat- 
ter is chosen so that it comprises information on special 
system-dependent behaviour (such as tendencies in the 
vicinity of system-specific singularities) to which the ex- 
ternal parameter is not sensitive. This method has been 
used successfully in Ref. [56] to calculate the transition 
line between the rung-singlet phase and a spontaneously 
dimerized phase for the spin ladder system including ring 
exchange. 

In this article we propose optimized perturbation the- 
ory (OPT), which is based on the principle of minimal 
sensitivity [49], as a particularly robust technique to si- 
multaneously extrapolate a large number of polynomials 
in an automatized fashion. 

For the general derivation of OPT we go back to the 
beginning of our perturbational approach where we as- 
sumed that the Hamiltonian can be split into an unper- 
turbed part U and a perturbation V 



H{x) = U + xV . 



(65) 



The fundamental idea of optimized perturbation theory 
(OPT) is to modify this splitting and to adjust this mod- 
ification by an additional control parameter a 



H{x; a) = (1 + a)U + xV - aU 

= (l + a)#(.i;a) 
H{x;d) = U + x{V + dU), X ''' 



(66) 



1 + a 



a 

X 



We consider x to be the new expansion parameter. The 
Hamiltonian H{x\a) = {1 + a)H{x\a) is identical to 
Hamiltonian H{x). Hence no energy depends on a in 
the exact result. Let us consider the gap A(a;) as generic 
example. It does not depend on a. But the truncated 
series expansion Atrunc(a;;a) resulting from H{x;a) will 
depend on a. We at least demand stationarity in this 
parameter. This is motivated by the idea of minimal 



sensitivity [49]. We write 



Atrunc(a;; a) = (1 + a)T \ A{x; d) , (67) 
where A(a;; d) denotes the energy resulting from H{x; d) 

n 

and T \ f{x) is the n**^ order Taylor expansion of /(x) 

X=Xo 

in x about x = xq- Requiring stationarity leads to the 
criterion 



9aAtruiic(a;;a)|a=aopt = 



(68) 



In general, Atrunc(a;; Oopt) converges faster than the cor- 
responding series expansions of A(.t), since the additional 
degree of freedom can be used to optimize the splitting 
into an unperturbed and a perturbing part [49]. In other 
words, the system has the freedom to choose the best 
splitting depending on the series under study. Moreover, 
in some cases a convergent series expansion can be en- 
forced by OPT even if the original series diverges. In 
Ref. [49], see also references therein, the harmonic os- 
cillator perturbed by a quartic potential is given as an 
example whose standard series expansion for the ground 
state energy diverges [63]. 

To be more specific, the series of A(i; a) is needed. We 
rewrite 



H{x;d) = U + x{V + dU) 



= {1 + dx) 



U + 



—V 



1 + dx 



(69) 



One clearly sees that the series of A in i can be obtained 
by expanding the expression 



A(i;a) = (1 -|-aS)A 



1 + dx 



(70) 



in X. Let \is assume that we had already calculated 
A(x) as a truncated series Atrunc(a;) from iJeflp- Then 
Atrunc(a;;a) is obtained by 



Atrunc(a;; a) = {l + a)T A{x; d) 

x=0 



(71) 



(1 + a)T 



x=0 



(1 -I- ai)Atrunc 



1+dx 



Atrunc(3^) before. Finally x and d are re-substitute by 
their definitions in Eq (66). 

In order to do all steps in one we introduce an auxiliary 
variable A for the derivation. Then, the Taylor expansion 
in X can be replaced by an expansion in A 



Xx 
1 + a 



and d 



a 

x 



(72) 
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where it is understood, that the final result is obtained 
for A = 1. This leads to 



z{x;a) 



(73) 



r 



(l + a(l- A))Atrunc 



A=0 



Xx 



1 + a(l - A) 



A=l 



The bottomline is that OPT can be used without com- 
puting any new coefficients. Only some straightforward 
computer algebra is needed. We take the direct expan- 
sions for the energy Atrunc(a;) obtained from the effective 
Hamiltonian and substitute and re-expand according to 
Eq. (73) to get the optimized expansions Atrunc(a^; ciopt) 
where Oopt is given by the minimal sensitivity crite- 
rion (68). 

From the discussion above it is clear that also other 
quantities A^runcix) like matrix elements of effective ob- 
servables can be optimized analogously 



At 



(74) 



T 



At 



A=0 



Ax 



l + a(l- A) 



J A=l 



The prefactor (1 + a(l — A)) is dropped since A does not 
depend on the global energy scale in H in contrast to 
A. Note that formula (73) can be used for any matrix 
element of the Hamiltonian. Formula (74) can be used 
for any dimensionless matrix element of an observable. 

The criterion of minimal sensitivity allows to elaborate 
further on the structure of flopt- We will show that 



■''opt 



(75) 



holds. Let i?trunc(a;; a) be the truncated series expansion 
of the quantity for which we want to find the optimum 
value flopt. In the following discussion R is an energy 
A or some matrix element A. To ease the notation we 
introduce the function 



gfy^ y\ _ |^'^trunc(w/w) for energies, 

' [Atrmciu/v) for matrix elements. 



(76) 



The derivative of g with respect to v is denoted by 
f{u,v) = dyg{u,v). The problem of calculating Oopt re- 
duces to 



0=daRtTunc{x; a) = T 



dag{Xx, 1 + a(l - A)) 



A=0 



T 



/(Aa;,l + a(l-A))(l-A) , 



A=0 



(77) 



where the notation for A = 1 in the end is suppressed for 
clarity. For the following argument it is important to see 
that 



T 



/(Aa;,l + a(l-A))(l-A) 



(78) 



A=0 



n-1 



fnX" + (1 - A)T 



/(Aa;,l + a(l-A)) , 



A=0 



holds, where /„ denotes the n^^ coefficient in the Taylor 
expansion of / with respect to A 



fn 



-idxrf{Xx,l+a{l-X)) . 



(79) 



For the final value A = 1, the second term on the right 
hand side of Eq. (78) vanishes. In addition, the structure 
in Eq. (79) is such, that with each derivation with respect 
to A we obtain either an a; or an a as internal derivative 
of the chain rule. Thus, setting A = 1 in the end, we find 



^a-^trunc('^) ^) — fn 



(80) 



to be a homogeneous polynomial in the variables x and 
a. In an n"^ order expansion the criterion of minimal 

sensitivity reads 



=a'opt 



i=0 



(81) 



which clearly shows, that we can always write aopt = 
Q!opta;. This proves the assertion (75). 

The proposed OPT procedure can be performed for all 
physical quantities of interest in particular for the matrix 
elements of Hi and H2. No new calculations are required. 
Instead, the plain series results can be promoted to OPT 
results by simple substitutions and re-expansion. 

In the application we modify the OPT idea slightly. 
We assume that there is an optimum splitting for each 
order, that means an optimum value for aopt, to do the 
perturbation expansion. This means that we do not 
adapt Qfopt to various different quantities, but we use 
one universal value depending only on the order of the 
series. This value is determined by simultaneously op- 
timizing some simple quantities like the one-triplon gap 
A(a;) = ui{x, k = tt) (Eq. (21)), the 5* = bound state en- 
ergy at K = TT {As=o{x)) and others with respect to the 
best Dlog-Pade approximant of these quantities. This ap- 
proach is based on the plausible assumption that all con- 
sidered quantities (here: energy levels) are governed by 
the same singularities reflecting the underlying physics. 
This idea is supported by the fact that all energy expan- 
sions we obtained start to deviate from their best extrap- 
olations at about the same value for x (cf. Fig. 14). Thus, 
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in contrast to the original spirit of the OPT method, we 
propose that aopt essentially depends on the model and 
the order of the expansions only, but not on the particu- 
lar quantity under study. 

In Fig. 14, A(x) = co{x,k = tt), Lo{x,k = 0) and 
As=o{x) are plotted as functions of x. The thin solid 
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FIG. 14. The elementary triplori gap A = aj(fe = tt), 
u>{k = 0) and the 5 = two-triplon gap as functions of 
the perturbation parameter x. For all energy levels the va- 
lidity of the plain scries results (thin lines) starts to break 
down at a; ~ 0.6. Various biased Dlog-Padc approximants 
(dpade in legend) are shown for each quantity as thick solid 
or long dashed lines. They yield the most reliable extrap- 
olations. The results obtained from optimized perturbation 
theory (OPT) with ciopt = 2.9a; arc depicted as thick short 
dashed lines. Except for u}{k = 0) there is no visible deviation 
from the Dlog-Pade results. For a; < 1 the optimized results 
can be used without appreciable loss of accuracy. 



lines correspond to the plain expansion results, reliable 
up to X « 0.6. Various Dlog-Pade approximants for each 
energy arc depicted by thick solid and thick long-dashed 
lines. They are biased by including the fact, that all en- 
ergies in the ladder system should grow linearly in x for 
a; — > oo. In that limit the system maps onto two de- 
coupled 5 = 1/2 chains whose energies E^, measured in 
units of the remaining coupling constant, are constants 
-Ei// -.^11= const.. Measuring these energies in units of Jj_, 
as we do, stipulates E^/J± = E^x/ => Eiy/J± ~ x for 
a; ^ 1 giving rise to the extrapolation bias. The thick 
short-dashed lines show the corresponding OPT-results 
with aopt = 2.9a;. The figure illustrates that it is possi- 
ble to choose a fixed aopt leading to a global improvement 
in all energy levels. For some levels the improvement is 
very good (e.g. As=o{x) or A (a;)), for others it is still 
reasonable good (e.g. lv{x,0)). 

The value for aopt depends on the order of the original 
truncated series. We found that aopt = 2.6a; gives best 



results for 13**^ order expansions. Thus, matrix elements 
of Hi are optimized with aopt = 2.9.T and those of H2 
with aopt = 2.6a;. 

It is a significant advantage that the OPT procedure as 
we use it is linear. Let Oa[-] denote the OPT procedure 
such that 



f{x; flopt) = Oa[f{x)] 



(82) 



is the optimized series obtained from the direct series 
f{x). Then, for a linear quantity F{x) = J^i^'-ifii^) one 
has 



F{x;aopt) = Oa 



^aifi{x) 

i 



(83) 



as long as all /j are given to the same order. For the two- 
particle interaction part H2 of i?off, for instance, this 
means that one can choose to optimize the matrix ele- 
ments {K,d'\H2\K,d) directly, or to optimize the two- 
particle interaction coefficients td-r,d' before the sums of 
the Fourier transform is carried out (see Eq. 34). This 
linearity is not ensured by Pade or Dlog-Pade extrapo- 
lations which represents a serious caveat in the practical 
use. 

Wc like to stress that OPT does not yield the best ap- 
proximants one can think of. It is rather a compromise 
between feasibility and quality. The OPT method rep- 
resents a very robust and smooth approximation scheme 
in the sense that none of the approximants diverges or 
produces unexpected pathologies. Its linearity makes it 
particular appropriate for the treatment of Fourier trans- 
formed matrix elements. 

Some of the matrix elements of Hi and H2 have been 
cross-checked with (Dlog-)Pade approximants leading to 
the conclusion that the proposed method is reliable up 
to X ~ 1 with a maximum error of about 5%. 

Probing the effect on the shape of the spectral densities 
by manually varying single matrix elements we find that 
the elements {Hi+H2)ij with i, j G {1,2, 3} influence the 
line shapes most. Naturally, matrix elements connecting 
short distances d have the biggest influence in a system 
with exponentially decreasing correlation lengths. Thus, 
to further improve our results, we replace these elements 
by the most reliable (Dlog-)Pade approximants of the 
underlying series expansions for each K considered. 



VI. SUMMARY 

In this article, the necessary details are given to under- 
stand how perturbative CUTs can be used to quantita- 
tively calculate the low- lying excitations of a certain class 
of many-particle systems. Particular emphasis is put 
on spectral densities of experimentally relevant observ- 
ables. Due to the finiteness of the perturbation order, the 
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method will work especially well for systems with short- 
range correlations. Furthermore, it must be possible to 
define suitable quasi-particles from which the whole spec- 
trum of the system under study can be constructed. The 
simplifications rendering high order results for the eigen 
energies possible arise from mapping the initial Hamil- 
tonian onto an effective one which conserves the number 
of particles. This enables separate calculations in the 0-, 
1-, 2-,... particle sectors. In each sector, only a few-body 
problem has to be solved. 

Effective observables representing measurement pro- 
cesses are obtained in a similar fashion. In general, they 
do not conserve the number of particles. But their ac- 
tion on relevant states can be classified according to the 
number of particles they inject into the system. 

The antiferromagnetic spin 1/2 ladder is used to illus- 
trate all steps of the calculations in detail. The pertur- 
bation is taken about the limit of isolated rungs. In this 
limit, rung-triplets are the elementary excitations, which 
suggests to name them "triplons" (cf. Rcf. [18]). They 
are suitable quasi-particles if the inter-rung interaction is 
switched on inducing a magnetic polarization cloud. 

The 0-, 1- and 2- particle sectors of the effective Hamil- 
tonian are studied separately. We address all possible 
difhculties including a discussion of the finite clusters 
needed to obtain the matrix elements for the infinite sys- 
tem by using to the linked cluster theorem. 

Four different observables are discussed for the ladder 
system relevant for neutron and light scattering experi- 
ments. We show in detail how the relevant quantities can 
be calculated to obtain the corresponding 1- and 2-triplon 
spectral densities for experiments at zero temperature. 

The perturbative CUT methods requires the extrapo- 
lation of a large number of quantities if the system is not 
very local. This is especially so for calculations in the 
two- and more-particles sectors. The article includes a 
general treatment of a novel extrapolation scheme (opti- 
mized perturbation theory, OPT) designed to simultane- 
ously extrapolate a large number of quantities. The OPT 
is introduced as a robust technique which does not nec- 
essarily render the best possible results. But it provides 
reliable results for the regimes of interest. The method 
is indispensable for situations where one needs automa- 
tized extrapolations, a task that can hardly be solved by 
standard techniques like Pade methods. 
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